function Ptf = CriterionOverallSubSectoral(Data,Ptf,Sectors)
        
    Quantile = Ptf.Quantiles;

    Intensity=Data.Intensity12;
    Scope=Data.Scope12;
    
    Intensity(isnan(Intensity)) = 0;

    Weights=Data.WeightsRescaled;
    tmp=Scope./Data.MktCap; 
    idx = ~isnan(tmp);
    Weights(isnan(tmp))=0;
    Weights(idx)=Weights(idx)/sum(Weights(idx));

        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    % For each subsector, determine nb of firms and sum of weights
        
    TableExcluded = table(Data.Name,Data.Sector,Data.PrimaryIndustry,Intensity,Weights,Data.x_OfCSO,'VariableNames',{'Name','Sector','PrimaryIndustry','Intensity','Weight','CSO'});
    [TableExcluded,indexTE] = sortrows(TableExcluded,4,'descend','MissingPlacement','last');

    ListSectors = unique(TableExcluded.Sector);
    nbS = size(ListSectors,1);

    NbFirmsbySectors = zeros(nbS,1);
    WeightExcluded = zeros(nbS,1);
    NbFirmsbySectorsAll = zeros(nbS,1);
    WeightAll = zeros(nbS,1);

    SectorData=[];
    
    for i=1:nbS

        idx = strcmp(TableExcluded.Sector,ListSectors(i));
        NbFirmsbySectors(i) = sum(idx);
        WeightExcluded(i) = sum(TableExcluded.Weight(idx));

        ListSubSectors = unique(TableExcluded.PrimaryIndustry(idx,:));
        nbSS = size(ListSubSectors,1);
        
        idx = strcmp(Data.Sector,ListSectors(i));
        NbFirmsbySectorsAll(i) = sum(idx);
        WeightAll(i) = sum(Data.WeightsRescaled(idx));

        for j=1:nbSS

            idxs = strcmp(TableExcluded.PrimaryIndustry,ListSubSectors(j));
            NbFirmsbySubSectorsAll(j) = sum(idxs);
            WeightAllSub(j) = sum(TableExcluded.Weight(idxs));

            SubSectorData = table(ListSectors(i),ListSubSectors(j),NbFirmsbySubSectorsAll(j),WeightAllSub(j));
                
            SectorData=[SectorData;SubSectorData];

        end
        
    end
    
    SectorData.Properties.VariableNames={'Sector','PrimaryIndustry','NbAll','WeightAll'};
    
    nbSSTot = size(SectorData,1);
    SectorData.NbExcluded = zeros(nbSSTot,1);
    SectorData.WeightExcluded = zeros(nbSSTot,1);    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    SumWeightExcluded=0;
    idxfirm=1;
    CriterionExcl = zeros(size(Intensity));
    SectorData.OK = ones(nbSSTot,1);
    
    while SumWeightExcluded < Quantile
    
        CandidateFirm = TableExcluded(idxfirm,:);
        
        CandidatePrimInd = CandidateFirm.PrimaryIndustry;
        
        idxSS = strcmp(SectorData.PrimaryIndustry,CandidatePrimInd);
        
        if (SectorData.WeightExcluded(idxSS) + CandidateFirm.Weight/SectorData.WeightAll(idxSS) > 0.5)  || (SectorData.OK(idxSS) == 0)
            SectorData.OK(idxSS) = 0;
            idxfirm=idxfirm+1;
            continue
        else
            SectorData.WeightExcluded(idxSS) = SectorData.WeightExcluded(idxSS) + CandidateFirm.Weight/SectorData.WeightAll(idxSS);
            SectorData.NbExcluded(idxSS) = SectorData.NbExcluded(idxSS) + 1;
            SumWeightExcluded = SumWeightExcluded + CandidateFirm.Weight;
            CriterionExcl(indexTE(idxfirm)) = 1;
            idxfirm=idxfirm+1;
        end
    
    end

    CriterionExcl = logical(CriterionExcl);
    WeightExcl = Weights(logical(CriterionExcl));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    % Reinvestment at the sectoral level
    
    nSS = size(SectorData,1);
    GICSS = Data.PrimaryIndustry;

    CriterionIncl = zeros(length(Intensity),1);
    pWeights = Weights;

    WeightIncluded=zeros(nSS,1);

    for i=1:nSS

       if SectorData.WeightExcluded(i) > 0
           
           idx= (categorical(GICSS)==SectorData{i,2});
           MktCap = Weights(idx)/sum(Weights(idx));
           tmp = [Intensity(idx)   MktCap   Weights(idx)];      
           tmp2 = sortrows(tmp,1,'ascend','MissingPlacement','last');    
           tmp2(:,2) = cumsum(tmp2(:,2))/sum(tmp2(:,2));
           indx = find(tmp2(:,2)<=0.5);
           if isempty(indx)
               indx=1;
           end
           ThresholdD = tmp2(indx(end),1);
           CriterionIncl(idx) = logical(Intensity(idx)<=ThresholdD);

           sumWeightExcl = SectorData.WeightExcluded(i) * SectorData.WeightAll(i);
           Weightsi =  Weights(idx);
           WeightIncl = Weightsi(logical(CriterionIncl(idx)));
           sumWeightIncl = sum(WeightIncl);

           Weightsi(logical(CriterionIncl(idx))) = WeightIncl * (1+sumWeightExcl/sumWeightIncl);
           pWeights(idx) = Weightsi;
       
           WeightIncluded(i) = sum(WeightIncl * sumWeightExcl/sumWeightIncl);

       end
       
    end
       
    keepallpWeights = pWeights;
    pWeights = pWeights(~CriterionExcl);
    
    CriterionExcl = logical(CriterionExcl);
    CriterionIncl = logical(CriterionIncl);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Ptf.NbFirmsExcluded = sum(CriterionExcl);
    Ptf.MktCapExcluded = sum(WeightExcl);
    Ptf.NbFirmsIncluded = sum(CriterionIncl);
    Ptf.MktCapIncluded = sum(WeightIncluded);
    
    keepScore = ~isnan(Scope);
    CleanWeights = Weights(keepScore);
    tmp=(Scope(keepScore))./Data.Revenues(keepScore);
    idx = ~isnan(tmp);
    InitWAIntensity = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx))*1000;
    disp('Initial WAIntensity')
    disp(InitWAIntensity)

    keepScore = ~isnan(Scope);
    CleanWeights = Weights(keepScore);
    tmp=(Scope(keepScore))./Data.MktCap(keepScore);
    idx = ~isnan(tmp);
    InitFootprint = sum(CleanWeights(idx) .* tmp(idx))*1000000 / sum(CleanWeights(idx));
    disp('Initial Footprint')
    disp(InitFootprint)

    disp('List of firms with highest carbon footprint')
    Texcl = table(Data.Name(CriterionExcl,:),Intensity(CriterionExcl,:),Weights(CriterionExcl,:),Scope(CriterionExcl,:),Data.MktCap(CriterionExcl,:),Data.x_OfCSO(CriterionExcl,:),Data.MarketValue(CriterionExcl,:),Data.Revenues(CriterionExcl,:),Data.Sector(CriterionExcl,:),Data.PrimaryIndustry(CriterionExcl,:),'VariableNames',{'Name','Intensity','Weight','Scope','MktCap','CSO','MarketValue','Revenues','Sector','PrimaryIndustry'});
    Texcls = sortrows(Texcl,2,'descend');
    tmp=Texcl.Scope./Texcl.Revenues;
    idx = ~isnan(tmp);
    ExcludedWAIntensity = sum(Texcl.Weight(idx) .* tmp(idx)) / sum(Texcl.Weight(idx));
    tmp=Texcl.Scope./Texcl.MktCap;
    idx = ~isnan(tmp);
    ExcludedFootprint = sum(Texcl.Weight(idx) .* tmp(idx)) / sum(Texcl.Weight(idx))*1000000;
    cumulex = table(cumsum(Texcls.Weight),'VariableNames',{'Cumul'});
    disp([Texcls cumulex])
    disp('Number of firms excluded')
    disp(size(Texcls,1))
    disp('Fraction of total equity excluded')
    disp(cumulex{end,:})
    disp('Excluded WAIntensity')
    disp(ExcludedWAIntensity*cumulex{end,:})
    disp('Excluded Footprint')
    disp(ExcludedFootprint*cumulex{end,:})
    
    disp('List of firms with lowest carbon footprint')
    Tincl = table(Data.Name(CriterionIncl,:),Intensity(CriterionIncl,:),keepallpWeights(CriterionIncl,:),Scope(CriterionIncl,:),Data.MktCap(CriterionIncl,:),Data.x_OfCSO(CriterionIncl,:),Data.MarketValue(CriterionIncl,:),Data.Revenues(CriterionIncl,:),Data.Sector(CriterionIncl,:),Weights(CriterionIncl,:),'VariableNames',{'Name','Intensity','Weight','Scope','MktCap','CSO','MarketValue','Revenues','Sector','OldWeight'});
    Tincls = sortrows(Tincl,2,'ascend');
    tmp=Tincl.Scope./Tincl.Revenues;
    idx = ~isnan(tmp);
    AddedWAIntensity = sum(Tincl.Weight(idx) .* tmp(idx)) / sum(Tincl.Weight(idx));
    tmp=Tincl.Scope./Tincl.MktCap;
    idx = ~isnan(tmp);
    AddedFootprint = sum(Tincl.Weight(idx) .* tmp(idx)) / sum(Tincl.Weight(idx))*1000000;
    cumulin = table(cumsum(Tincls.Weight),'VariableNames',{'Cumul'});
    disp([Tincls cumulin])
    disp('Number of firms included')
    disp(size(Tincls,1))
    disp('Fraction of total equity included')
    disp(sum(WeightIncluded))
    disp('Added WAIntensity')
    disp(AddedWAIntensity*cumulin{end,:})
    disp('Added Footprint')
    disp(AddedFootprint*cumulin{end,:})

    DataAfterExclusion = Data;
    DataAfterExclusion(CriterionExcl,:) = [];
    NewWeights = pWeights / sum(pWeights);

    keepScore = ~isnan(DataAfterExclusion.Scope12);
    CleanWeights = NewWeights(keepScore) / sum(NewWeights(keepScore));
    tmp=(DataAfterExclusion.Scope12(keepScore))./DataAfterExclusion.Revenues(keepScore);
    idx = ~isnan(tmp);
    Ptf.WAIntensity = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx));
    disp('Final WAIntensity')
    disp(Ptf.WAIntensity)

    tmp=(DataAfterExclusion.Scope12(keepScore))./DataAfterExclusion.MktCap(keepScore);
    idx = ~isnan(tmp);
    Ptf.Footprint = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx))*1000000;
    disp('Final Footprint')
    disp(Ptf.Footprint)
    
    tmp1=DataAfterExclusion.Scope12(keepScore) .* DataAfterExclusion.x_OfCSO(keepScore);
    tmp2=DataAfterExclusion.Revenues(keepScore) .* DataAfterExclusion.x_OfCSO(keepScore);
    idx = ~isnan(tmp1.*tmp2);
    Ptf.Intensity = sum(tmp1(idx)) / sum(tmp2(idx))/1;

    tmp=DataAfterExclusion.TRIp1(keepScore)./DataAfterExclusion.TRI(keepScore)-1;
    idx = ~isnan(tmp);
    Ptf.AnnualReturn = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx));
